rm(list = ls(all = TRUE))
setwd("C:\\Users\\lenovo\\Desktop\\R")
#从线性模型到广义线性模型的推广
#7.1 二分类Logistic回归
#7.1.1 Logistic回归模型
p <- seq(from = 0, to = 1, by = .01)
odds <- p/(1-p)
plot(log(odds), p, type = "l", col = "blue", 
     ylab = "Probability", las =1)
abline(h = 0.5, lty = "dashed")
abline(v = 0, lty = "dashed")
#用于估计给定暴露水平时结果事件发生的概率

#7.1.2 Logistic回归示例
#1 数据准备
library(MASS)
data(birthwt)
str(birthwt)
#查看ptl（早产次数）、ftv（怀孕初期探访医生的次数）
library(epiDisplay)
tab1(birthwt$ptl)
#无早产史的就占了84%，转换为二因子更为合理
tab1(birthwt$ftv)
#将探访医生次数超过一次的孕妇合并成一类，转换为3因子

library(dplyr)
birthweight <- birthwt %>%
  mutate(race = factor(race, labels = c("white", "black", "other")),
         smoke = factor(smoke, labels = c("no", "yes")),
         ptl = ifelse(ptl > 0, "1+", ptl),
         ptl = factor(ptl),
         ht = factor(ht, labels = c("no", "yes")),
         ui = factor(ui, labels = c("no", "yes")),
         ftv = ifelse(ftv > 1, "2+", ftv),
         ftv = factor(ftv))
str(birthweight)

#2 模型的建立
#binomial代表Logistic回归
glm1 <- glm(low~age + lwt + race + smoke + ptl + ht + ui + ftv,
            family = binomial,
            data = birthweight)
summary(glm1)
#两个离差平方和的差值为134.67-195.48 = 39.19；计算p值
#Number of Fisher Scoring iterations: 4 
#这个指标一般不管；但如果大于25，说明模型过于复杂
pchisq(39.19, df = 188-178, lower.tail = FALSE)
#也可以用anova得到；但是anova一般用于比较包含不同自变量的模型

#3 自变量的筛选
drop1(glm1)
glm2 <- step(glm1)
summary(glm2)

#4 模型的比较

#glm2剔除了age与ftv
#对于两个嵌套模型（glm2的自变量为glm1的子集），
#使用似然比检验（Likelihood Ratio Test，LRT）
#LLR = 2（lnL1-lnL2）
#L1为复杂模型（含较多自变量的模型）的似然值
#L2为简单模型的似然值
lnL1 <- as.numeric(logLik(glm1))
lnL2 <- as.numeric(logLik(glm2))
#因为变量age包含一个参数，
#变量ftv包含两个参数（三个水平的因子需要两个哑变量）
#模型参数个数差为3，自由度只差为df=3
pchisq(2*(lnL1-lnL2), df = 3, lower.tail = FALSE)
anova(glm1, glm2, test = "Chisq")
anova(glm1, glm2, test = "LRT")

#也可以使用AIC比较两个或多个模型
AIC(glm1, glm2)

#5 回归系数的解释
#先查看回归系数
coef(glm2)
#把结果指数转换为优势比
exp(coef(glm2))
#lwt约为0.98，如果lwt增加10（lb），优势增加0.98^10

#获取回归系数的置信区间
confint(glm2)
#获取优势比的置信区间
exp(confint(glm2))

#6 预测
#虚构一个案例
newdata <- data.frame(lwt = 120, race = "black",
                      smoke = "yes", ptl = "0", 
                      ht = "yes", ui = "no")
logit <- predict(glm2, newdata = newdata)
logit
#转换为概率
exp(logit)/(1+exp(logit))
#其分娩低体重儿的概率约为88%

#或者直接用predict函数
predict(glm2, newdata = newdata, type = "response")

#7 模型的检查

#Hosmer-Lemeshow检验可以用于判断观测值与预测值之间差异的显著性
#install.packages("ResourceSelection")
library(ResourceSelection)
hoslem.test(birthweight$low, fitted(glm2))
#不能拒绝关于模型很好地拟合了数据的假设

#8 模型结果的汇总输出
library(epiDisplay)
logistic.display(glm2)

write.csv(logistic.display(glm2)$table, file = "model.csv")


#7.1.3 表格数据的Logistic回归
#数据输入,注意输入的顺序
dat.array <- array(c(136,57,107,151,63,44,63,265),
                   dim = c(2,2,2),
                   dimnames = list(smoke = c("no","yes"),
                                   drink = c("no","yes"),
                                   outcome = c("control","case")))
dat.table <- as.table(dat.array)
dat.table
dat <- as.data.frame(dat.table)
dat

#模型构建
#weights设为频数
mod1 <- glm(outcome~smoke+drink,family = binomial,
            weights = Freq, data = dat)
summary(mod1)

#或者，把结果变量outcome分为两列
library(tidyr)
dat.wide <- pivot_wider(dat, names_from = outcome, values_from = Freq)
dat.wide
#模型构建
mod2 <- glm(cbind(case,control)~smoke + drink,
            family = binomial, data = dat.wide)
summary(mod2)

#转换成原始数据后建立模型
dat.original <- dat[rep(1:nrow(dat),dat$Freq),-4]
nrow(dat.original)
mod3 <- glm(outcome ~ smoke+ drink, family = binomial,
            data = dat.original)
summary(mod3)

#7.2 条件Logistic回归——配对资料

#条件Logistic回归模型的不同之处在常数项
#不同匹配组的常数项可以不同
#而每个变量在不同匹配组的相应相同——即回归系数相同

#载入数据
library(epiDisplay)
data(VC1to1)
str(VC1to1)
head(VC1to1)

#建立模型
library(survival)
clogit1 <- clogit(case~smoking+alcohol+rubber+strata(matset),
                  data = VC1to1)
clogit1
#无统计学意义

drop1(clogit1, test = "Chisq")
#变量smoking对应的p值最大，先剔除
clogit2 <- clogit(case~alcohol+rubber+strata(matset),
                  data = VC1to1)
clogit2
drop1(clogit2, test = "Chisq")
#rubber仍然有较大的p值

clogit3 <- clogit(case~alcohol+strata(matset),
                  data = VC1to1)
clogit3

#比较
AIC(clogit1,clogit2,clogit3)
summary(clogit3)

#汇总结果、输出表格
clogistic.display(clogit3)
write.csv(clogistic.display(clogit3)$table, file = "clogit.csv")

#7.3 无序多分类Logistic回归——基线分类logit模型

#加载数据
library(epiDisplay)
data(Ectopic)
str(Ectopic)

#默认把分类变量outc的第一个水平EP（“宫外孕患者”）作为参照组
#可以换（似乎失败了）
#Ectopic$outcIA <- relevel(Ectopic$outc, ref="IA")

#install.packages("nnet")
library(nnet)
multi1 <- multinom(outc~hia, data = Ectopic)
summary(multi1)
#该研究为一个病例对照研究，常数项没有意义
#关注变量hia的系数

#由回归系数与标准误获得每个系数的z值
st <- summary(multi1)$standard.errors
z <- coef(multi1)/st
z

#获取p值
p.values <- pnorm(abs(z), lower.tail = FALSE)*2
p.values

#回归系数的95%置信区间
confint(multi1)

#优势比
exp(coef(multi1))
#优势比的95%置信区间
exp(confint(multi1))

#较简单地获取以上数据
library(epiDisplay)
mlogit.display(multi1)

#把变量gravi加入模型
multi2 <- multinom(outc~hia + gravi, data = Ectopic)
mlogit.display(multi2)

#模型比较,离差减小
lrtest(multi1,multi2)
AIC(multi1,multi2)

#?Ectopic

#7.4 有序Logistic回归

#录入数据
dat <- array(c(10,7,19,6,0,2,7,5,1,5,6,16),
             dim = c(2,2,3),
             dimnames = list(method = c("old","new"),
                             sex = c("male","female"),
                             outcome = c("effectless","effective","recover")))

dat <- as.table(dat)
data1 <- as.data.frame(dat)
data1
str(data1)

#把结果变量转换为有序因子
#由于函数默认以字母顺序从低到高排列因子的各个水平
#此次无需指定顺序
data1$outcome <- ordered(data1$outcome)
data1$outcome
#对比无序变量
data1$method

#建立累积优势Logistic模型
library(MASS)
polr1 <- polr(outcome~sex+method, weights = Freq, data = data1)
summary(polr1)

#也可以把数据框data1转换成逐条记录的数据后建立模型
data1
data2 <- data1[rep(1:nrow(data1),data1$Freq),-4]
data2
head(data2)
str(data2)
polr2 <- polr(outcome~sex+method, data = data2)
summary(polr2)

#α1=1.8128，α2=2.6672，β1=1.319，β2=1.797
#ln（P1/P2+P3）=1.8128-1.319X1-1.797X2
#ln（P1+P2/P3）=2.6672-1.319X1-1.797X2

#回归系数指数化
exp(coef(polr1))
exp(confint(polr1))

#本例中的OR是“无效对有效和痊愈”或“无效和有效对痊愈”
#两个变量的OR值都大于1,且置信区间都不包含1
#女性患者的疗效显著优于男性患者 sexfemale 1.352945 11.33812
#新疗法的疗效显著优于旧疗效 methodnew 2.458614 15.77240


#筛选自变量
#step
#drop1

#对模型的平行性假设做检验
#install.packages("VGAM")
library(VGAM)
polr.p <- vglm(outcome~sex + method,
               family = cumulative(parallel = TRUE),
               data = data2)
polr.np <- vglm(outcome~sex+method,
                family = cumulative(parallel = FALSE),
                data = data2)
#避免与epiDisplay中的lrtest重名
VGAM::lrtest(polr.p, polr.np)
#两模型差异没有统计学意义，因此可认为模型满足平行性假设条件

#导出
library(epiDisplay)
ordinal.or.display(polr1)
write.csv(ordinal.or.display(polr1),file = "polr.csv")


